{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Wall time: 280 ms\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "from osgeo import gdal\n",
    "import numpy as np\n",
    "import os \n",
    "from glob import glob\n",
    "from math import ceil\n",
    "import time"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 边界\n",
    "def GetExtent(infile):\n",
    "    ds = gdal.Open(infile)\n",
    "    geotrans = ds.GetGeoTransform()\n",
    "    xsize = ds.RasterXSize\n",
    "    ysize = ds.RasterYSize\n",
    "    min_x,max_y = geotrans[0],geotrans[3]\n",
    "    max_x,min_y = geotrans[0]+xsize*geotrans[1],geotrans[3]+ysize*geotrans[5]\n",
    "    ds = None\n",
    "    return min_x,max_y,max_x,min_y"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def compress(path, target_path,method=\"LZW\"): #\n",
    "    \"\"\"使用gdal进行文件压缩，\n",
    "          LZW方法属于无损压缩，\n",
    "          效果非常给力，4G大小的数据压缩后只有三十多M\"\"\"\n",
    "    dataset = gdal.Open(path)\n",
    "    driver = gdal.GetDriverByName('GTiff')\n",
    "    driver.CreateCopy(target_path, dataset, strict=1, options=[\"TILED=YES\", \"COMPRESS={0}\".format(method)])\n",
    "    del dataset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def RasterMosaic(file_list,outpath):\n",
    "    Open = gdal.Open\n",
    "    min_x,max_y,max_x,min_y=GetExtent(file_list[0])\n",
    "    for infile in file_list:\n",
    "        minx,maxy,maxx,miny = GetExtent(infile)\n",
    "        min_x,min_y = min(min_x,minx),min(min_y,miny)\n",
    "        max_x,max_y = max(max_x,maxx),max(max_y,maxy)\n",
    "    \n",
    "    in_ds = Open(file_list[0])\n",
    "    in_band=in_ds.GetRasterBand(1)\n",
    "    geotrans = list(in_ds.GetGeoTransform())\n",
    "    width,height = geotrans[1],geotrans[5]\n",
    "    columns = ceil((max_x-min_x)/width)#列数\n",
    "    rows = ceil((max_y-min_y)/(-height))#行数\n",
    "    \n",
    "    outfile = outpath+file_list[0][:2]+'.tif' #结果文件名，可自行修改\n",
    "    driver=gdal.GetDriverByName('GTiff')\n",
    "    # out_ds=driver.Create(outfile,columns,rows,1,in_band.DataType)\n",
    "    out_ds=driver.Create(outfile,columns,rows,1,1)\n",
    "    out_ds.SetProjection(in_ds.GetProjection())\n",
    "    geotrans[0]=min_x#更正左上角坐标\n",
    "    geotrans[3]=max_y\n",
    "    out_ds.SetGeoTransform(geotrans)\n",
    "    out_band=out_ds.GetRasterBand(1)\n",
    "    inv_geotrans=gdal.InvGeoTransform(geotrans)\n",
    "\n",
    "    for in_fn in file_list:\n",
    "        in_ds=Open(in_fn)\n",
    "        in_gt=in_ds.GetGeoTransform()\n",
    "        offset=gdal.ApplyGeoTransform(inv_geotrans,in_gt[0],in_gt[3])\n",
    "        x,y=map(int,offset)\n",
    "\n",
    "        data=in_ds.GetRasterBand(1).ReadAsArray()\n",
    "\n",
    "        data = (data / 93 * 256).astype(np.uint8)\n",
    "\n",
    "        out_band.WriteArray(data,x,y)#x，y是开始写入时左上角像元行列号\n",
    "    del in_ds,out_band,out_ds\n",
    "    return outfile"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# path = r'J:/backup' #该文件夹下存放了待拼接的栅格\n",
    "path = r'E:/datasets/spacenet-dataset/spacenet/SN6_buildings/train/AOI_11_Rotterdam/PS-RGB' #该文件夹下存放了待拼接的栅格\n",
    "os.chdir(path)\n",
    "raster_list = sorted(glob('*.tif')) #读取文件夹下所有tif数据\n",
    "result = RasterMosaic(raster_list,outpath = r'E:/datasets/spacenet-dataset/spacenet/SN6_buildings/train/AOI_11_Rotterdam/PS-RGB' ) #拼接栅格\n",
    "# compress(result,target_path = r'E:/datasets/unname/test_public/AOI_11_Rotterdam/SAR-Intensity/Global_.tif') #压缩栅格"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "path = r'J:/backup' #该文件夹下存放了待拼接的栅格\n",
    "os.chdir(path)\n",
    "raster_list = sorted(glob('*.tif')) #读取文件夹下所有tif数据\n",
    "result = RasterMosaic(raster_list,outpath = r'J:/backup/Global' ) #拼接栅格\n",
    "compress(result,target_path = r'J:/backup/Global.tif') #压缩栅格"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3.6.15 ('SOMatch36')",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.15"
  },
  "orig_nbformat": 4,
  "vscode": {
   "interpreter": {
    "hash": "30d4882d9e1aba50c89ababc08c8aefb10a5319dbc4507c8f7697d7a0def6c41"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
